function Ptf = CriterionOverallSectoral(Data,Ptf,Sectors)
        
    Quantile = Ptf.Quantiles;

    Intensity=Data.Intensity12;
    Scope=Data.Scope12;
    
    Intensity(isnan(Intensity)) = 0;

    Weights=Data.WeightsRescaled;
    tmp=Scope./Data.MktCap; 
    idx = ~isnan(tmp);
    Weights(isnan(tmp))=0;
    Weights(idx)=Weights(idx)/sum(Weights(idx));

        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    % Exclusion at the portfolio level
    
    tmp = [Intensity Weights];
    tmp2 = sortrows(tmp,1,'descend','MissingPlacement','last');
    tmp2(:,2) = cumsum(tmp2(:,2));
    indx = find(tmp2(:,2)<=Quantile);
    if isempty(indx)
        indx=1;
    end
    ThresholdU = tmp2(indx(end),1);

    CriterionExcl = logical(Intensity>=ThresholdU);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Identify sectors

    ExcludedData = Data(CriterionExcl,:);

    ListSectors = unique(ExcludedData.Sector);
    nbS = size(ListSectors,1);

    NbFirmsbySectors = zeros(nbS,1);
    WeightExcluded = zeros(nbS,1);
    NbFirmsbySectorsAll = zeros(nbS,1);
    WeightAll = zeros(nbS,1);

    for i=1:nbS

        idx = strcmp(ExcludedData.Sector,ListSectors(i));
        NbFirmsbySectors(i) = sum(idx);
        WeightExcluded(i) = sum(ExcludedData.WeightsRescaled(idx));

        idx = strcmp(Data.Sector,ListSectors(i));
        NbFirmsbySectorsAll(i) = sum(idx);
        WeightAll(i) = sum(Data.WeightsRescaled(idx));

    end

    TSectors = table(ListSectors,NbFirmsbySectors,WeightExcluded);
    TSectors = sortrows(TSectors,3,{'descend'});
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    % Reinvestment at the sectoral level
    
    nS = size(TSectors,1);
    GICSS = Data.Sector;

    CriterionIncl = zeros(length(Intensity),1);
    pWeights = Weights;
    WeightIncluded = zeros(nS,1);

    for i=1:nS

       sumWeightExcl = TSectors{i,3};
        
       idx= (categorical(GICSS)==TSectors{i,1});
       MktCap = Weights(idx)/sum(Weights(idx));
       tmp = [Intensity(idx) MktCap];

       WeightExcludedWithinSector = sumWeightExcl/sum(Weights(idx));
       if WeightExcludedWithinSector < 0.5
           Quantilei = WeightExcludedWithinSector; % I overweight the same fraction of MC as I exclude firms
       else
           Quantilei = 1-WeightExcludedWithinSector; % I overweight the same fraction of MC as I exclude firms
       end
       
       tmp2 = sortrows(tmp,1,'ascend','MissingPlacement','last');    
       tmp2(:,2) = cumsum(tmp2(:,2));
       indx = find(tmp2(:,2)<=Quantilei);
       if isempty(indx)
           indx=1;
       end
       ThresholdD = tmp2(indx(end),1);

       CriterionIncl(idx) = logical(Intensity(idx)<=ThresholdD);
       
       CriterionIncl2 = zeros(length(Intensity),1);
       CriterionIncl2(idx) = logical(Intensity(idx)<=ThresholdD);

       idx2 = CriterionIncl2.*(Intensity<=ThresholdD);

       Weightsi = Weights(idx);
       WeightIncl = Weightsi(logical(CriterionIncl(idx)));
       sumWeightIncl = sum(WeightIncl);
       WeightIncluded(i) = sumWeightIncl;

       Weightsi(logical(CriterionIncl(idx))) = WeightIncl * (1+sumWeightExcl/sumWeightIncl);

       pWeights(idx) = Weightsi;
       
       Data.WeightsAdded(logical(idx2)) = WeightIncl * (sumWeightExcl/sumWeightIncl);
                  
    end

    keepallpWeights = pWeights;
    keepallpWeights(CriterionExcl) = 0;
    pWeights = pWeights(~CriterionExcl);
    
    CriterionExcl = logical(CriterionExcl);
    CriterionIncl = logical(CriterionIncl);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Ptf.NbFirmsExcluded = sum(CriterionExcl);
    Ptf.MktCapExcluded = sum(WeightExcluded);
    Ptf.NbFirmsIncluded = sum(CriterionIncl);
    Ptf.MktCapIncluded = sum(WeightIncluded);
    
    keepScore = ~isnan(Scope);
    CleanWeights = Weights(keepScore);
    tmp=(Scope(keepScore))./Data.Revenues(keepScore);
    idx = ~isnan(tmp);
    InitWAIntensity = sum(CleanWeights(idx) .* tmp(idx)) / sum(CleanWeights(idx))*1000;
    disp('Initial WAIntensity')
    disp(InitWAIntensity)

    keepScore = ~isnan(Scope);
    CleanWeights = Weights(keepScore);
    tmp=(Scope(keepScore))./Data.MktCap(keepScore);
    idx = ~isnan(tmp);
    InitFootprint = sum(CleanWeights(idx) .* tmp(idx))*1000000 / sum(CleanWeights(idx));
    disp('Initial Footprint')
    disp(InitFootprint)

    disp('List of firms with highest carbon footprint')
    Texcl = table(Data.Name(CriterionExcl,:),Intensity(CriterionExcl,:),Weights(CriterionExcl,:),Scope(CriterionExcl,:),Data.MktCap(CriterionExcl,:),Data.x_OfCSO(CriterionExcl,:),Data.MarketValue(CriterionExcl,:),Data.Revenues(CriterionExcl,:),Data.Sector(CriterionExcl,:),'VariableNames',{'Name','Intensity','Weight','Scope','MktCap','CSO','MarketValue','Revenues','Sector'});
    Texcls = sortrows(Texcl,2,'descend');
    tmp=Texcl.Scope./Texcl.Revenues;
    idx = ~isnan(tmp);
    ExcludedWAIntensity = sum(Texcl.Weight(idx) .* tmp(idx)) / sum(Texcl.Weight(idx));
    tmp=Texcl.Scope./Texcl.MktCap;
    idx = ~isnan(tmp);
    ExcludedFootprint = sum(Texcl.Weight(idx) .* tmp(idx)) / sum(Texcl.Weight(idx))*1000000;
    cumulex = table(cumsum(Texcls.Weight),'VariableNames',{'Cumul'});
    disp([Texcls cumulex])
    disp('Number of firms excluded')
    disp(size(Texcls,1))
    disp('Exclusion threshold')
    disp(ThresholdU)
    disp('Fraction of total equity excluded')
    disp(cumulex{end,:})
    disp('Excluded WAIntensity')
    disp(ExcludedWAIntensity*cumulex{end,:})
    disp('Excluded Footprint')
    disp(ExcludedFootprint*cumulex{end,:})
    
    disp('List of firms with lowest carbon footprint')
    Tincl = table(Data.Name(CriterionIncl,:),Intensity(CriterionIncl,:),Data.WeightsAdded(CriterionIncl,:),Scope(CriterionIncl,:),Data.MktCap(CriterionIncl,:),Data.x_OfCSO(CriterionIncl,:),Data.MarketValue(CriterionIncl,:),Data.Revenues(CriterionIncl,:),Data.Sector(CriterionIncl,:),Data.WeightsRescaled(CriterionIncl,:),keepallpWeights(CriterionIncl,:),'VariableNames',{'Name','Intensity','Weight','Scope','MktCap','CSO','MarketValue','Revenues','Sector','Initial_Weight','Final_Weight'});
    Tincls = sortrows(Tincl,2,'ascend');
    tmp=Tincl.Scope./Tincl.Revenues;
    idx = ~isnan(tmp);
    AddedWAIntensity = sum(Tincl.Weight(idx) .* tmp(idx)) / sum(Tincl.Weight(idx));
    tmp=Tincl.Scope./Tincl.MktCap;
    idx = ~isnan(tmp);
    AddedFootprint = sum(Tincl.Weight(idx) .* tmp(idx)) / sum(Tincl.Weight(idx))*1000000;
    cumulin = table(cumsum(Tincls.Weight),'VariableNames',{'Cumul'});
    disp([Tincls cumulin])
    disp('Number of firms included')
    disp(size(Tincls,1))
    disp('Overweighting threshold')
    disp(ThresholdD)
    disp('Fraction of total equity included')
    disp(cumulin{end,:})
    disp('Added WAIntensity')
    disp(AddedWAIntensity*cumulin{end,:})
    disp('Added Footprint')
    disp(AddedFootprint*cumulin{end,:})

    DataAfterExclusion = Data;
    DataAfterExclusion(CriterionExcl,:) = [];
    NewWeights = pWeights / sum(pWeights);

    keepScore = ~isnan(DataAfterExclusion.Scope12);
    CleanWeights = NewWeights(keepScore) / sum(NewWeights(keepScore));
    tmp=(DataAfterExclusion.Scope12(keepScore))./DataAfterExclusion.Revenues(keepScore);
    idx = ~isnan(tmp);
    Ptf.WAIntensity = sum(CleanWeights(idx) .* tmp(idx)) / sum(CleanWeights(idx));
    disp('Final WAIntensity')
    disp(Ptf.WAIntensity)
   
    tmp=(DataAfterExclusion.Scope12(keepScore))./DataAfterExclusion.MktCap(keepScore);
    idx = ~isnan(tmp);
    Ptf.Footprint = sum(CleanWeights(idx) .* tmp(idx)) / sum(CleanWeights(idx))*1000000;
    disp('Final Footprint')
    disp(Ptf.Footprint)

    tmp1=DataAfterExclusion.Scope12(keepScore) .* DataAfterExclusion.x_OfCSO(keepScore);
    tmp2=DataAfterExclusion.Revenues(keepScore) .* DataAfterExclusion.x_OfCSO(keepScore);
    idx = ~isnan(tmp1.*tmp2);
    Ptf.Intensity = sum(tmp1(idx)) / sum(tmp2(idx))/1;

    tmp=DataAfterExclusion.TRIp1(keepScore)./DataAfterExclusion.TRI(keepScore)-1;
    idx = ~isnan(tmp);
    Ptf.AnnualReturn = sum(CleanWeights(idx) .* tmp(idx)) / sum(CleanWeights(idx));
        
   